home *** CD-ROM | disk | FTP | other *** search
/ Aminet 6 / Aminet 6 - June 1995.iso / Aminet / gfx / 3d / irit50src.lha / irit5 / symb_lib / symbsply.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-09-23  |  43.9 KB  |  1,173 lines
  1. /******************************************************************************
  2. * SymbSPly.c - Adaptive surface to polygons approximation.              *
  3. *******************************************************************************
  4. * Written by Gershon Elber, March. 93.                          *
  5. ******************************************************************************/
  6.  
  7. #include "symb_loc.h"
  8. #include "geomat3d.h"
  9.  
  10. typedef enum {
  11.     UMIN_BOUNDARY = 1,
  12.     UMAX_BOUNDARY = 2,
  13.     VMIN_BOUNDARY = 4,
  14.     VMAX_BOUNDARY = 8
  15. } PatchBoundaryType;
  16.  
  17. typedef struct PointDependStruct {
  18.     struct PointDependStruct *Pnext;
  19.     struct RectangularStruct *Rect;
  20.     int OtherIndex, Index;
  21. } PointDependStruct;
  22.  
  23. typedef struct RectangularStruct {
  24.     struct RectangularStruct *Pnext;
  25.     CagdPType Point[4];    /* 4 corners, not always on Srf due to "black holes".*/
  26.     CagdVType Normal[4];
  27.     CagdUVType UV[4];
  28.     PointDependStruct *Depend[4];  /* Points that depends on these 4 points. */
  29. } RectangularStruct;
  30.  
  31. typedef struct EdgeStruct {
  32.     struct EdgeStruct *Pnext;
  33.     struct RectangularStruct *Rect;
  34.     int Index;
  35. } EdgeStruct;
  36.  
  37. typedef struct {
  38.     struct EdgeStruct *Left;        /* UMin */
  39.     struct EdgeStruct *Right;        /* UMax */
  40.     struct EdgeStruct *Top;        /* VMin */
  41.     struct EdgeStruct *Bottom;        /* VMax */
  42. } ConnectivityStruct;
  43.  
  44. static int
  45.     GlblFoundPartialShare = FALSE;
  46. static RectangularStruct
  47.     *GlblRectList = NULL;
  48. static CagdSrfStruct
  49.     *GlblCrvtrSqrSrf = NULL,
  50.     *GlblUIsoCrvtrSqrSrf = NULL,
  51.     *GlblVIsoCrvtrSqrSrf = NULL,
  52.     *GlblOrigSrf = NULL,
  53.     *GlblNormalSrf = NULL;
  54.  
  55. static CagdRType SymbMaxDistCrvLine(CagdCrvStruct *Crv);
  56. static CagdRType ExtremumNodeValue(CagdSrfStruct *Srf,
  57.                    CagdRType UMin, CagdRType UMax,
  58.                    CagdRType VMin, CagdRType VMax,
  59.                    CagdSrfDirType Dir,
  60.                    CagdRType *MaxVal);
  61. static ConnectivityStruct SymbSrf2OptimalPolygonsAux(CagdSrfStruct *Srf,
  62.                 int SubdivDepth,
  63.                 CagdRType Tolerance,
  64.                 SymbPlSubdivStrategyType SubdivDirStrategy,
  65.                 SymbPlErrorFuncType SrfPolyApproxError);
  66. static ConnectivityStruct MergeConnectivityEdge(
  67.                 ConnectivityStruct *Connectivity1,
  68.                 ConnectivityStruct *Connectivity2,
  69.                 CagdSrfDirType Dir);
  70. static void UpdateBlackHoles(EdgeStruct *L1,
  71.                  EdgeStruct *L2,
  72.                  CagdSrfDirType Dir);
  73. static int UpdateOneBlackHole(RectangularStruct *Rect1,
  74.                   int Index1a,
  75.                   RectangularStruct *Rect2,
  76.                   int Index2a,
  77.                   CagdSrfDirType Dir);
  78. static void MakePointDepend(RectangularStruct *DependRect,
  79.                 int DependIndex,
  80.                 RectangularStruct *MasterRect,
  81.                 int Master1Index,
  82.                 int Master2Index);
  83. static void PropagateDependency(RectangularStruct *Rect, int Index);
  84. static void ProjectPointOnLine(CagdPType Pt,
  85.                    CagdPType Pt1,
  86.                    CagdPType Pt2,
  87.                    CagdPType Nrml,
  88.                    CagdPType Nrml1,
  89.                    CagdPType Nrml2);
  90. static EdgeStruct *ChainEdgeLists(EdgeStruct *L1, EdgeStruct *L2);
  91. static CagdPolygonStruct *MakeTriangle(CagdPType Pt1,
  92.                        CagdVType Nrml1,
  93.                        CagdUVType UV1,
  94.                        CagdPType Pt2,
  95.                        CagdVType Nrml2,
  96.                        CagdUVType UV2,
  97.                        CagdPType Pt3,
  98.                        CagdVType Nrml3,
  99.                        CagdUVType UV3,
  100.                        CagdBType ComputeUV,
  101.                        CagdBType ComputeNrmls);
  102. static EdgeStruct *MakeEdge(RectangularStruct *Rect, int Index);
  103.  
  104. /*****************************************************************************
  105. * DESCRIPTION:                                                               *
  106. * Routine to bound the maximal deviation of a curve from a line connecting   *
  107. * its end points. Works for closed curves as well.                 *
  108. *                                                                            *
  109. * PARAMETERS:                                                                *
  110. *   Crv:        To compute a bound on its straightness.                      *
  111. *                                                                            *
  112. * RETURN VALUE:                                                              *
  113. *   CagdRType:    Straightness deviation measure.                            *
  114. *****************************************************************************/
  115. static CagdRType SymbMaxDistCrvLine(CagdCrvStruct *Crv)
  116. {
  117.     int i;
  118.     CagdRType *R,
  119.     MaxVal = 0.0;
  120.     CagdPtStruct Pt1, Pt2;
  121.     CagdCrvStruct *Line, *Diff, *DiffSqr;
  122.  
  123.     CagdCoerceToE3(Pt1.Pt, Crv -> Points, 0, Crv -> PType);
  124.     CagdCoerceToE3(Pt2.Pt, Crv -> Points, Crv -> Length - 1, Crv -> PType);
  125.  
  126.     Line = CagdMergePtPt(&Pt1, &Pt2);
  127.     if (CAGD_IS_BSPLINE_CRV(Crv)) {
  128.     CagdCrvStruct
  129.         *CTmp = CnvrtBezier2BsplineCrv(Line);
  130.  
  131.     CagdCrvFree(Line);
  132.     Line = CTmp;
  133.     BspKnotAffineTrans(Line -> KnotVector, Line -> Length + Line -> Order,
  134.                Crv -> KnotVector[0] - Line -> KnotVector[0],
  135.                (Crv -> KnotVector[Crv -> Length +
  136.                           Crv -> Order - 1] -
  137.                 Crv -> KnotVector[0]) /
  138.                (Line -> KnotVector[Line -> Length +
  139.                            Line -> Order - 1] -
  140.                 Line -> KnotVector[0]));
  141.     }
  142.     Diff = SymbCrvSub(Crv, Line);
  143.     CagdCrvFree(Line);
  144.     DiffSqr = SymbCrvDotProd(Diff, Diff);
  145.     CagdCrvFree(Diff);
  146.     Diff = CagdCoerceCrvTo(DiffSqr, CAGD_PT_E1_TYPE);
  147.     CagdCrvFree(DiffSqr);
  148.     for (i = 0, R = Diff -> Points[1]; i < Diff -> Length; i++, R++) {
  149.     if (MaxVal < *R)
  150.         MaxVal = *R;
  151.     }
  152.     return sqrt(MaxVal);
  153. }
  154.  
  155. /*****************************************************************************
  156. * DESCRIPTION:                                                               M
  157. * Routine to compute the scalar field of k1^2 + k2^2 (k1, k2 are principal   M
  158. * curvatures) for the surface Srf, into GlblCrvtrSqrSrf.             M
  159. *   This scalar field is used by SymbSrf2OptPolysCurvatureError function.    M
  160. *                                                                            *
  161. * PARAMETERS:                                                                M
  162. *   Srf:        To compute the curvature bound for as an optional preprocess M
  163. *        for function SymbSrf2OptPolysCurvatureError.             M
  164. *                                                                            *
  165. * RETURN VALUE:                                                              M
  166. *   void                                                                     M
  167. *                                                                            *
  168. * KEYWORDS:                                                                  M
  169. *   SymbSrf2OptPolysCurvatureErrorPrep, polygonization, surface             M
  170. *   approximation                                 M
  171. *****************************************************************************/
  172. void SymbSrf2OptPolysCurvatureErrorPrep(CagdSrfStruct *Srf)
  173. {
  174.     if (CAGD_IS_BEZIER_SRF(Srf))
  175.     GlblOrigSrf = CnvrtBezier2BsplineSrf(Srf);
  176.     else
  177.     GlblOrigSrf = CagdSrfCopy(Srf);
  178.  
  179.     GlblCrvtrSqrSrf = SymbSrfCurvatureUpperBound(Srf);
  180. }
  181.  
  182. /*****************************************************************************
  183. * DESCRIPTION:                                                               M
  184. * Routine to estimate the curvature of the patch using k1^2 + k2^2.          M
  185. *   Assumes the availability of the GlblCrvtrSqrSrf for Srf.             M
  186. *   This estimate is too loose and in fact is not recommended!             M
  187. *                                                                            *
  188. * PARAMETERS:                                                                M
  189. *   Srf:         To estimate curvature for.                                  M
  190. *   Dir:         Currently not used.                                         M
  191. *                                                                            *
  192. * RETURN VALUE:                                                              M
  193. *   CagdRType:    Curvature estimated.                                       M
  194. *                                                                            *
  195. * KEYWORDS:                                                                  M
  196. *   SymbSrf2OptPolysCurvatureError, polygonization, surface approximation    M
  197. *****************************************************************************/
  198. CagdRType SymbSrf2OptPolysCurvatureError(CagdSrfStruct *Srf,
  199.                      CagdSrfDirType Dir)
  200. {
  201.     RealType UMin, UMax, VMin, VMax, *R, MaxCurvatureSqr, Size, RetVal, TmpR;
  202.     CagdSrfStruct *TSrf1, *TSrf2;
  203.     CagdCrvStruct *AuxCrv, *Crv;
  204.     CagdBBoxStruct BBox;
  205.     int i,
  206.     Neighbors = AttrGetIntAttrib(Srf -> Attr, "_Neighbors");
  207.  
  208.     if ((RetVal = AttrGetRealAttrib(Srf -> Attr, "_SrfCurvature")) <
  209.                             IP_ATTR_BAD_REAL)
  210.     return RetVal;
  211.  
  212.     CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
  213.  
  214.     TSrf1 = CagdSrfRegionFromSrf(GlblCrvtrSqrSrf, UMin, UMax, CAGD_CONST_U_DIR);
  215.     TSrf2 = CagdSrfRegionFromSrf(TSrf1, VMin, VMax, CAGD_CONST_V_DIR);
  216.     CagdSrfFree(TSrf1);
  217.     TSrf1 = CagdCoerceSrfTo(TSrf2, CAGD_PT_E1_TYPE);
  218.     CagdSrfFree(TSrf2);
  219.  
  220.     for (i = 0, R = TSrf1 -> Points[1], MaxCurvatureSqr = 0.0;
  221.      i < TSrf1 -> ULength * TSrf1 -> VLength;
  222.      i++, R++)
  223.     if (MaxCurvatureSqr < *R)
  224.         MaxCurvatureSqr = *R;
  225.     CagdSrfFree(TSrf1);
  226.  
  227.     CagdSrfBBox(Srf, &BBox);
  228.     for (i = 0, Size = 0.0; i < 3; i++)
  229.     Size += BBox.Max[i] - BBox.Min[i];
  230.     
  231.     RetVal = sqrt(MaxCurvatureSqr) * Size;
  232.     
  233.     /* Now compute bounds on boundary, if on boundary. */
  234.     if (Neighbors & UMIN_BOUNDARY) {
  235.     AuxCrv = CagdCrvFromMesh(Srf, 0, CAGD_CONST_U_DIR),
  236.     Crv = CagdCrvRegionFromCrv(AuxCrv, VMin, VMax);
  237.     TmpR = SymbMaxDistCrvLine(Crv);
  238.     RetVal = MAX(RetVal, TmpR);
  239.     CagdCrvFree(AuxCrv);
  240.     CagdCrvFree(Crv);
  241.     }
  242.     if (Neighbors & UMAX_BOUNDARY) {
  243.     AuxCrv = CagdCrvFromMesh(Srf, Srf -> ULength - 1,
  244.                  CAGD_CONST_U_DIR);
  245.     Crv = CagdCrvRegionFromCrv(AuxCrv, VMin, VMax);
  246.     TmpR = SymbMaxDistCrvLine(Crv);
  247.     RetVal = MAX(RetVal, TmpR);
  248.     CagdCrvFree(AuxCrv);
  249.     CagdCrvFree(Crv);
  250.     }
  251.     if (Neighbors & VMIN_BOUNDARY) {
  252.     AuxCrv = CagdCrvFromMesh(Srf, 0, CAGD_CONST_V_DIR);
  253.     Crv = CagdCrvRegionFromCrv(AuxCrv, UMin, UMax);
  254.     TmpR = SymbMaxDistCrvLine(Crv);
  255.     RetVal = MAX(RetVal, TmpR);
  256.     CagdCrvFree(AuxCrv);
  257.     CagdCrvFree(Crv);
  258.     }
  259.     if (Neighbors & VMAX_BOUNDARY) {
  260.     AuxCrv = CagdCrvFromMesh(Srf, Srf -> VLength - 1,
  261.                  CAGD_CONST_V_DIR);
  262.     Crv = CagdCrvRegionFromCrv(AuxCrv, UMin, UMax);
  263.     TmpR = SymbMaxDistCrvLine(Crv);
  264.     RetVal = MAX(RetVal, TmpR);
  265.     CagdCrvFree(AuxCrv);
  266.     CagdCrvFree(Crv);
  267.     }
  268.  
  269.     AttrSetRealAttrib(&Srf -> Attr, "_SrfCurvature", RetVal);
  270.     return RetVal;
  271. }
  272.  
  273. /*****************************************************************************
  274. * DESCRIPTION:                                                               M
  275. * Routine to estimate the curvature of the patch using a bilinear approx.    M
  276. *                                                                            *
  277. * PARAMETERS:                                                                M
  278. *   Srf:         To estimate curvature for.                                  M
  279. *   Dir:         Currently not used.                                         M
  280. *                                                                            *
  281. * RETURN VALUE:                                                              M
  282. *   CagdRType:    Curvature estimated.                                       M
  283. *                                                                            *
  284. * KEYWORDS:                                                                  M
  285. *   SymbSrf2OptPolysBilinPolyError, polygonization, surface approximation    M
  286. *****************************************************************************/
  287. CagdRType SymbSrf2OptPolysBilinPolyError(CagdSrfStruct *Srf,
  288.                      CagdSrfDirType Dir)
  289. {
  290.     int i;
  291.     RealType UMin, UMax, VMin, VMax, *R, MaxDiffSqr, RetVal, Bilin2PolyErr;
  292.     CagdPtStruct Pt00, Pt01, Pt10, Pt11;
  293.     CagdSrfStruct *BilinSrf, *TSrf2,
  294.     *TSrf1 = CagdSrfCopy(Srf);
  295.     PlaneType Plane;
  296.  
  297.     CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
  298.     R = CagdSrfEval(Srf, UMin, VMin);
  299.     CagdCoerceToE3(Pt00.Pt, &R, -1, Srf -> PType);
  300.     R = CagdSrfEval(Srf, UMin, VMax);
  301.     CagdCoerceToE3(Pt01.Pt, &R, -1, Srf -> PType);
  302.     R = CagdSrfEval(Srf, UMax, VMin);
  303.     CagdCoerceToE3(Pt10.Pt, &R, -1, Srf -> PType);
  304.     R = CagdSrfEval(Srf, UMax, VMax);
  305.     CagdCoerceToE3(Pt11.Pt, &R, -1, Srf -> PType);
  306.  
  307.     if (CGPlaneFrom3Points(Plane, Pt00.Pt, Pt01.Pt, Pt10.Pt))
  308.     Bilin2PolyErr = CGDistPointPlane(Pt11.Pt, Plane) / 2;
  309.     else
  310.     Bilin2PolyErr = 0.0;
  311.  
  312.     BilinSrf = CagdBilinearSrf(&Pt00, &Pt10, &Pt01, &Pt11);
  313.  
  314.     CagdMakeSrfsCompatible(&TSrf1, &BilinSrf, TRUE, TRUE, TRUE, TRUE);
  315.  
  316.     TSrf2 = SymbSrfSub(TSrf1, BilinSrf);
  317.     CagdSrfFree(TSrf1);
  318.     CagdSrfFree(BilinSrf);
  319.     TSrf1 = SymbSrfDotProd(TSrf2, TSrf2);
  320.     CagdSrfFree(TSrf2);
  321.  
  322.     TSrf2 = CagdCoerceSrfTo(TSrf1, CAGD_PT_E1_TYPE);
  323.     CagdSrfFree(TSrf1);
  324.  
  325.     for (i = 0, R = TSrf2 -> Points[1], MaxDiffSqr = 0.0;
  326.      i < TSrf2 -> ULength * TSrf2 -> VLength;
  327.      i++, R++)
  328.     if (MaxDiffSqr < *R)
  329.         MaxDiffSqr = *R;
  330.     CagdSrfFree(TSrf2);
  331.  
  332.     RetVal = sqrt(MaxDiffSqr) + Bilin2PolyErr;
  333.     
  334.     AttrSetRealAttrib(&Srf -> Attr, "_SrfCurvature", RetVal);
  335.     return RetVal;
  336. }
  337.  
  338. /*****************************************************************************
  339. * DESCRIPTION:                                                               M
  340. * Routine to precompute the scalar field of kn^u and kn^v (the normal        M
  341. * curvatures in the iso parametric directions).                     M
  342. *   These scalar fields are used to determined the prefered subdivision         M
  343. * location of Srf.                                 M
  344. *                                                                            *
  345. * PARAMETERS:                                                                M
  346. *   Srf:      To compute the curvature bound in the isoparametric direction. M
  347. *                                                                            *
  348. * RETURN VALUE:                                                              M
  349. *   void                                                                     M
  350. *                                                                            *
  351. * KEYWORDS:                                                                  M
  352. *   SymbSrf2OptPolysIsoDirCurvatureErrorPrep, polygonization, surface         M
  353. *   approximation                                 M
  354. *****************************************************************************/
  355. void SymbSrf2OptPolysIsoDirCurvatureErrorPrep(CagdSrfStruct *Srf)
  356. {
  357.     if (CAGD_IS_BEZIER_SRF(Srf))
  358.     GlblOrigSrf = CnvrtBezier2BsplineSrf(Srf);
  359.     else
  360.     GlblOrigSrf = CagdSrfCopy(Srf);
  361.  
  362.     GlblUIsoCrvtrSqrSrf = SymbSrfIsoDirNormalCurvatureBound(Srf,
  363.                                 CAGD_CONST_U_DIR);
  364.     GlblVIsoCrvtrSqrSrf = SymbSrfIsoDirNormalCurvatureBound(Srf,
  365.                                 CAGD_CONST_V_DIR);
  366. }
  367.  
  368. /*****************************************************************************
  369. * DESCRIPTION:                                                               M
  370. * Routine to convert a single surface to a set of triangles approximating    M
  371. * it.                                         M
  372. *   FineNess is controlled via a function that returns an error measure      M
  373. * SrfPolyApproxError that is guaranteed to be less than Tolerance.         M
  374. *                                                                            *
  375. * PARAMETERS:                                                                M
  376. *   Srf:               To convert and approximate using triangles.           M
  377. *   Tolerance:         Accuracy control.                                     M
  378. *   SubdivDirStrategy: Alternatively in U and V, direction that minimizes    M
  379. *               the error, etc.                         M
  380. *   SrfPolyApproxErr:  Using bilinear curvature estimate, k1^2 + k2^2         M
  381. *               estimate, etc.                         M
  382. *   ComputeNormals:    Do we want normals to be computed as well?            M
  383. *   FourPerFlat:       If TRUE, four triangle per flat surface patch are     M
  384. *                      created, otherwise only two.                 M
  385. *   ComputeUV:         Do we want UV parameter values with the vertices of   M
  386. *                      the triangles?                         M
  387. *                                                                            *
  388. * RETURN VALUE:                                                              M
  389. *   CagdPolygonStruct *:  Resulting polygons that approximates Srf.          M
  390. *                                                                            *
  391. * KEYWORDS:                                                                  M
  392. *   SymbSrf2OptimalPolygons, approximation, conversion                       M
  393. *****************************************************************************/
  394. CagdPolygonStruct *SymbSrf2OptimalPolygons(CagdSrfStruct *Srf,
  395.                 CagdRType Tolerance,
  396.                 SymbPlSubdivStrategyType SubdivDirStrategy,
  397.                 SymbPlErrorFuncType SrfPolyApproxErr,
  398.                 CagdBType ComputeNormals,
  399.                 CagdBType FourPerFlat,
  400.                 CagdBType ComputeUV)
  401. {
  402.     int i,
  403.     NewSrf = FALSE;
  404.     CagdPolygonStruct
  405.     *Polys = NULL;
  406.  
  407.     GlblNormalSrf = ComputeNormals ? SymbSrfNormalSrf(Srf) : NULL;
  408.  
  409.     if (CAGD_IS_BEZIER_SRF(Srf)) {
  410.     NewSrf = TRUE;
  411.     Srf = CnvrtBezier2BsplineSrf(Srf);
  412.     }
  413.  
  414.     GlblRectList = NULL;
  415.     GlblFoundPartialShare = FALSE;
  416.  
  417.     AttrSetIntAttrib(&Srf -> Attr, "_Neighbors",
  418.              UMIN_BOUNDARY | UMAX_BOUNDARY |
  419.              VMIN_BOUNDARY | VMAX_BOUNDARY);
  420.  
  421.     SymbSrf2OptimalPolygonsAux(Srf, 0, Tolerance, SubdivDirStrategy,
  422.                    SrfPolyApproxErr);
  423.  
  424.     /* Convert every rectangular domain into two or four polygons. */
  425.     while (GlblRectList != NULL) {
  426.     RectangularStruct
  427.         *Rect = GlblRectList;
  428.  
  429.     if (FourPerFlat) {
  430.         CagdPType MiddlePoint;
  431.         CagdVType MiddleNormal;
  432.         CagdUVType MiddleUV;
  433.         CagdRType *R;
  434.  
  435.         MiddleUV[0] = (Rect -> UV[0][0] + Rect -> UV[2][0]) / 2.0;
  436.         MiddleUV[1] = (Rect -> UV[0][1] + Rect -> UV[2][1]) / 2.0;
  437.         R = CagdSrfEval(Srf, MiddleUV[0], MiddleUV[1]);
  438.         CagdCoerceToE3(MiddlePoint, &R, -1, Srf -> PType);
  439.         CagdEvaluateSurfaceVecField(MiddleNormal, GlblNormalSrf,
  440.                     MiddleUV[0], MiddleUV[1]);
  441.  
  442.         for (i = 0; i < 4; i++) {
  443.         CagdPolygonStruct
  444.             *Poly = MakeTriangle(Rect -> Point[i],
  445.                          Rect -> Normal[i],
  446.                          Rect -> UV[i],
  447.                      Rect -> Point[(i + 1) % 4],
  448.                          Rect -> Normal[(i + 1) % 4],
  449.                          Rect -> UV[(i + 1) % 4],
  450.                      MiddlePoint,
  451.                          MiddleNormal,
  452.                          MiddleUV,
  453.                      ComputeUV, ComputeNormals);
  454.  
  455.         if (Poly != NULL)
  456.             LIST_PUSH(Poly, Polys);
  457.         }
  458.     }
  459.     else {
  460.         CagdPolygonStruct
  461.         *Poly1 = MakeTriangle(Rect -> Point[0],
  462.                           Rect -> Normal[0],
  463.                           Rect -> UV[0],
  464.                       Rect -> Point[1],
  465.                           Rect -> Normal[1],
  466.                           Rect -> UV[1],
  467.                       Rect -> Point[2],
  468.                           Rect -> Normal[2],
  469.                           Rect -> UV[2],
  470.                       ComputeUV, ComputeNormals),
  471.         *Poly2 = MakeTriangle(Rect -> Point[0],
  472.                           Rect -> Normal[0],
  473.                           Rect -> UV[0],
  474.                       Rect -> Point[2],
  475.                           Rect -> Normal[2],
  476.                           Rect -> UV[2],
  477.                       Rect -> Point[3],
  478.                           Rect -> Normal[3],
  479.                           Rect -> UV[3],
  480.                       ComputeUV, ComputeNormals);
  481.  
  482.         if (Poly1 != NULL)
  483.         LIST_PUSH(Poly1, Polys);
  484.         if (Poly2 != NULL)
  485.         LIST_PUSH(Poly2, Polys);
  486.     }
  487.  
  488.     GlblRectList = GlblRectList -> Pnext;
  489.  
  490.     /* Free the rectangular structure. */
  491.     for (i = 0; i < 4; i++) {
  492.         PointDependStruct
  493.         *Depend = Rect -> Depend[i];
  494.  
  495.         while (Depend != NULL) {
  496.         PointDependStruct
  497.             *NextDepend = Depend -> Pnext;
  498.  
  499.         IritFree((VoidPtr) Depend);
  500.         Depend = NextDepend;
  501.         }
  502.     }
  503.     IritFree((VoidPtr) Rect);
  504.     }
  505.  
  506.     CagdSrfFree(GlblNormalSrf);
  507.     if (NewSrf)
  508.     CagdSrfFree(Srf);
  509.  
  510.     if (GlblCrvtrSqrSrf != NULL) {
  511.     CagdSrfFree(GlblCrvtrSqrSrf);
  512.     GlblCrvtrSqrSrf = NULL;
  513.     }
  514.     if (GlblUIsoCrvtrSqrSrf != NULL) {
  515.     CagdSrfFree(GlblUIsoCrvtrSqrSrf);
  516.     GlblUIsoCrvtrSqrSrf = NULL;
  517.     }
  518.     if (GlblVIsoCrvtrSqrSrf != NULL) {
  519.     CagdSrfFree(GlblVIsoCrvtrSqrSrf);
  520.     GlblVIsoCrvtrSqrSrf = NULL;
  521.     }
  522.     if (GlblOrigSrf != NULL) {
  523.     CagdSrfFree(GlblOrigSrf);
  524.     GlblOrigSrf = NULL;
  525.     }
  526.  
  527.     return Polys;
  528. }
  529.  
  530. /*****************************************************************************
  531. * DESCRIPTION:                                                               *
  532. * Computs the extremum curvature location to subdivide at.             *
  533. *                                                                            *
  534. * PARAMETERS:                                                                *
  535. *   Srf:             Surface to consider.                         *
  536. *   UMin, UMax, VMin, VMax:     Domain to consider in Srf.                   *
  537. *   Dir:                        Direction of subdivision planned.            *
  538. *   MaxVal:                     Where extremum curvature value is saved.     *
  539. *                                                                            *
  540. * RETURN VALUE:                                                              *
  541. *   CagdRType:    The subdivision parameter in direction Dir suggested       *
  542. *****************************************************************************/
  543. static CagdRType ExtremumNodeValue(CagdSrfStruct *Srf,
  544.                    CagdRType UMin,
  545.                    CagdRType UMax,
  546.                    CagdRType VMin,
  547.                    CagdRType VMax,
  548.                    CagdSrfDirType Dir,
  549.                    CagdRType *MaxVal)
  550. {
  551.     CagdRType *UV;
  552.     CagdSrfStruct
  553.     *TSrf1 = CagdSrfRegionFromSrf(Srf, UMin, UMax, CAGD_CONST_U_DIR),
  554.     *TSrf2 = CagdSrfRegionFromSrf(TSrf1, VMin, VMax, CAGD_CONST_V_DIR);
  555.  
  556.     CagdSrfFree(TSrf1);
  557.     TSrf1 = CagdCoerceSrfTo(TSrf2, CAGD_PT_E1_TYPE);
  558.     CagdSrfFree(TSrf2);
  559.  
  560.     UV = BspSrfMaxCoefParam(TSrf1, 1, MaxVal);
  561.     CagdSrfFree(TSrf1);
  562.  
  563.     if (APX_EQ(UV[0], UMin))
  564.     UV[0] = (2.0 * UMin + UMax) / 3.0;
  565.     if (APX_EQ(UV[0], UMax))
  566.     UV[0] = (2.0 * UMax + UMin) / 3.0;
  567.     if (APX_EQ(UV[1], VMin))
  568.     UV[1] = (2.0 * VMin + VMax) / 3.0;
  569.     if (APX_EQ(UV[1], VMax))
  570.     UV[1] = (2.0 * VMax + VMin) / 3.0;
  571.  
  572.     return Dir == CAGD_CONST_U_DIR ? UV[0] : UV[1];
  573. }
  574.  
  575. /*****************************************************************************
  576. * DESCRIPTION:                                                               *
  577. * Routine to convert a single surface to set of rectangular domains         *
  578. * approximating it.                                 *
  579. *   FineNess is control via a function that returns an error measure         *
  580. * SrfPolyApproxError that is guaranteed to be less than Tolerance.         *
  581. *   The returned Connectivity structure is used to fill in "black holes".    *
  582. *                                                                            *
  583. * PARAMETERS:                                                                *
  584. *   Srf:               To convert and approximate using triangles.           *
  585. *   SubdivDepth:       Depth of recursion (subdivision) of this function.    *
  586. *   Tolerance:         Accuracy control.                                     *
  587. *   SubdivDirStrategy: Alternatively in U and V, direction that minimizes    *
  588. *               the error, etc.                         *
  589. *   SrfPolyApproxErr:  Using bilinear curvature estimate, k1^2 + k2^2         *
  590. *               estimate, etc.                         *
  591. *                                                                            *
  592. * RETURN VALUE:                                                              *
  593. *   ConnectivityStruct:   Subdivided regions with adjacency topology.        *
  594. *****************************************************************************/
  595. static ConnectivityStruct SymbSrf2OptimalPolygonsAux(CagdSrfStruct *Srf,
  596.                 int SubdivDepth,
  597.                 CagdRType Tolerance,
  598.                 SymbPlSubdivStrategyType SubdivDirStrategy,
  599.                 SymbPlErrorFuncType SrfPolyApproxErr)
  600. {
  601.     CagdRType UMin, UMax, VMin, VMax, u, v;
  602.     int UOrder = Srf -> UOrder,
  603.     VOrder = Srf -> VOrder,
  604.     ULength = Srf -> ULength,
  605.     VLength = Srf -> VLength,
  606.     Neighbors = AttrGetIntAttrib(Srf -> Attr, "_Neighbors");
  607.     CagdBType
  608.     HasUDiscont = Srf -> UKnotVector != NULL &&
  609.             BspKnotC1Discont(Srf -> UKnotVector, UOrder, ULength, &u),
  610.     HasVDiscont = Srf -> VKnotVector != NULL &&
  611.             BspKnotC1Discont(Srf -> VKnotVector, VOrder, VLength, &v);
  612.  
  613.     CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
  614.  
  615. #ifdef DEBUG_DOMAIN
  616.     fprintf(stderr,
  617.         "Umin = %lf, Umax = %lf, Vmin = %lf, Vmax = %lf (%d, %d)\n",
  618.         UMin, UMax, VMin, VMax, HasUDiscont, HasVDiscont);
  619. #endif /* DEBUG_DOMAIN */
  620.  
  621.     if (!(HasUDiscont || HasVDiscont) &&
  622.     (UMax - UMin < EPSILON * 5 ||
  623.      VMax - VMin < EPSILON * 5 ||
  624.      SrfPolyApproxErr(Srf, CAGD_NO_DIR) < Tolerance)) {
  625.     int i;
  626.     ConnectivityStruct Connectivity;
  627.     RectangularStruct
  628.         *Rect = (RectangularStruct *)
  629.         IritMalloc(sizeof(RectangularStruct));
  630.  
  631.     CagdSrfDomain(Srf,
  632.               &Rect -> UV[0][0], &Rect -> UV[2][0],
  633.               &Rect -> UV[0][1], &Rect -> UV[2][1]);
  634.     Rect -> UV[1][0] = Rect -> UV[2][0];
  635.     Rect -> UV[1][1] = Rect -> UV[0][1];
  636.     Rect -> UV[3][0] = Rect -> UV[0][0];
  637.     Rect -> UV[3][1] = Rect -> UV[2][1];
  638.  
  639.     for (i = 0; i < 4; i++) {
  640.         CagdRType
  641.         *R = CagdSrfEval(Srf, Rect -> UV[i][0], Rect -> UV[i][1]);
  642.  
  643.         CagdCoerceToE3(Rect -> Point[i], &R, -1, Srf -> PType);
  644.  
  645.         CagdEvaluateSurfaceVecField(Rect -> Normal[i], GlblNormalSrf,
  646.                     Rect -> UV[i][0], Rect -> UV[i][1]);
  647.  
  648.         Rect -> Depend[i] = NULL;
  649.     }
  650.     LIST_PUSH(Rect, GlblRectList);
  651.  
  652.     Connectivity.Left = MakeEdge(Rect, 3);
  653.     Connectivity.Right = MakeEdge(Rect, 1);
  654.     Connectivity.Top = MakeEdge(Rect, 0);
  655.     Connectivity.Bottom = MakeEdge(Rect, 2);
  656.     return Connectivity;
  657.     }
  658.     else {
  659.     int NeighborU1 = Neighbors & ~UMAX_BOUNDARY,
  660.         NeighborU2 = Neighbors & ~UMIN_BOUNDARY,
  661.         NeighborV1 = Neighbors & ~VMAX_BOUNDARY,
  662.         NeighborV2 = Neighbors & ~VMIN_BOUNDARY;
  663.     CagdRType UDiv, VDiv, ErrU1, ErrU2, ErrV1, ErrV2,
  664.         MaxUIsoCrvtr, MaxVIsoCrvtr;
  665.     CagdSrfStruct *SrfU1, *SrfU2, *SrfV1, *SrfV2;
  666.     ConnectivityStruct Connectivity, Connectivity1, Connectivity2;
  667.     CagdSrfDirType SubdivDir;
  668.  
  669.     if (HasUDiscont) {
  670.         NeighborU1 = Neighbors | UMAX_BOUNDARY;
  671.         NeighborU2 = Neighbors | UMIN_BOUNDARY;
  672.         UDiv = u;
  673.     }
  674.     else if (GlblUIsoCrvtrSqrSrf != NULL)
  675.         UDiv = ExtremumNodeValue(GlblVIsoCrvtrSqrSrf,
  676.                      UMin, UMax, VMin, VMax,
  677.                      CAGD_CONST_U_DIR, &MaxVIsoCrvtr);
  678.     else if (ULength > UOrder)
  679.         UDiv = Srf -> UKnotVector[(Srf -> ULength + Srf -> UOrder) / 2];
  680.     else
  681.         UDiv = (UMin + UMax) / 2;
  682.  
  683.     if (HasVDiscont) {
  684.         NeighborV1 = Neighbors | VMAX_BOUNDARY;
  685.         NeighborV2 = Neighbors | VMIN_BOUNDARY;
  686.         VDiv = v;
  687.     }
  688.     else if (GlblVIsoCrvtrSqrSrf != NULL)
  689.         VDiv = ExtremumNodeValue(GlblUIsoCrvtrSqrSrf,
  690.                      UMin, UMax, VMin, VMax,
  691.                      CAGD_CONST_V_DIR, &MaxUIsoCrvtr);
  692.     else if (VLength > VOrder)
  693.         VDiv = Srf -> VKnotVector[(Srf -> VLength + Srf -> VOrder) / 2];
  694.     else
  695.         VDiv = (VMin + VMax) / 2;
  696.     
  697.     SrfU1 = CagdSrfSubdivAtParam(Srf, UDiv, CAGD_CONST_U_DIR);
  698.     SrfU2 = SrfU1 -> Pnext;
  699.     AttrSetIntAttrib(&SrfU1 -> Attr, "_Neighbors", NeighborU1);
  700.     AttrSetIntAttrib(&SrfU2 -> Attr, "_Neighbors", NeighborU2);
  701.     
  702.     SrfV1 = CagdSrfSubdivAtParam(Srf, VDiv, CAGD_CONST_V_DIR);
  703.     SrfV2 = SrfV1 -> Pnext;
  704.     AttrSetIntAttrib(&SrfV1 -> Attr, "_Neighbors", NeighborV1);
  705.     AttrSetIntAttrib(&SrfV2 -> Attr, "_Neighbors", NeighborV2);
  706.  
  707.     if (HasUDiscont)
  708.         SubdivDir = CAGD_CONST_U_DIR;
  709.     else if (HasVDiscont)
  710.         SubdivDir = CAGD_CONST_V_DIR;
  711.     else {
  712.         if (GlblUIsoCrvtrSqrSrf != NULL && GlblVIsoCrvtrSqrSrf != NULL)
  713.             SubdivDir = MaxUIsoCrvtr > MaxVIsoCrvtr ? CAGD_CONST_V_DIR
  714.                             : CAGD_CONST_U_DIR;
  715.  
  716.         switch (SubdivDirStrategy) {
  717.         case SYMB_SUBDIV_STRAT_MIN_MAX:
  718.             ErrU1 = SrfPolyApproxErr(SrfU1, CAGD_CONST_U_DIR);
  719.             ErrU2 = SrfPolyApproxErr(SrfU2, CAGD_CONST_U_DIR);
  720.             ErrV1 = SrfPolyApproxErr(SrfV1, CAGD_CONST_V_DIR);
  721.             ErrV2 = SrfPolyApproxErr(SrfV2, CAGD_CONST_V_DIR);
  722.             SubdivDir = MAX(ErrU1, ErrU2) > MAX(ErrV1, ErrV2) ?
  723.             CAGD_CONST_V_DIR : CAGD_CONST_U_DIR;
  724.             break;
  725.         default:
  726.         case SYMB_SUBDIV_STRAT_ALTERNATE:
  727.             SubdivDir = SubdivDepth % 2 ?
  728.             CAGD_CONST_U_DIR : CAGD_CONST_V_DIR;
  729.             break;
  730.         }
  731.     }
  732.  
  733.     if (SubdivDir == CAGD_CONST_V_DIR) {
  734.         /* Use the V subdivision. */
  735.         Connectivity1 = SymbSrf2OptimalPolygonsAux(SrfV1, SubdivDepth + 1,
  736.                                Tolerance,
  737.                                SubdivDirStrategy,
  738.                                SrfPolyApproxErr);
  739.         Connectivity2 = SymbSrf2OptimalPolygonsAux(SrfV2, SubdivDepth + 1,
  740.                                Tolerance,
  741.                                SubdivDirStrategy,
  742.                                SrfPolyApproxErr);
  743.  
  744.         Connectivity = MergeConnectivityEdge(&Connectivity1,
  745.                          &Connectivity2,
  746.                          CAGD_CONST_V_DIR);
  747.     }
  748.     else {
  749.         /* Use the U subdivision. */
  750.         Connectivity1 = SymbSrf2OptimalPolygonsAux(SrfU1, SubdivDepth + 1,
  751.                                Tolerance,
  752.                                SubdivDirStrategy,
  753.                                SrfPolyApproxErr);
  754.         Connectivity2 = SymbSrf2OptimalPolygonsAux(SrfU2, SubdivDepth + 1,
  755.                                Tolerance,
  756.                                SubdivDirStrategy,
  757.                                SrfPolyApproxErr);
  758.  
  759.         Connectivity = MergeConnectivityEdge(&Connectivity1,
  760.                          &Connectivity2,
  761.                          CAGD_CONST_U_DIR);
  762.     }
  763.     CagdSrfFree(SrfU1);
  764.     CagdSrfFree(SrfU2);
  765.     CagdSrfFree(SrfV1);
  766.     CagdSrfFree(SrfV2);
  767.     return Connectivity;
  768.     }
  769. }
  770.  
  771. /*****************************************************************************
  772. * DESCRIPTION:                                                               *
  773. * Merges two connectivity information into one, destructively.             *
  774. *                                                                            *
  775. * PARAMETERS:                                                                *
  776. *   Connectivity1, Connectivity2:  The two regions to merge.             *
  777. *   Dir:                           Direction of merge. Either U or V.         *
  778. *                                                                            *
  779. * RETURN VALUE:                                                              *
  780. *   ConnectivityStruct:   Merged structure.                     *
  781. *****************************************************************************/
  782. static ConnectivityStruct MergeConnectivityEdge(
  783.                     ConnectivityStruct *Connectivity1,
  784.                     ConnectivityStruct *Connectivity2,
  785.                     CagdSrfDirType Dir)
  786. {
  787.     ConnectivityStruct Connectivity;
  788.  
  789.     switch (Dir) {
  790.     case CAGD_CONST_U_DIR:
  791.         Connectivity.Left = Connectivity1 -> Left;
  792.         Connectivity.Right = Connectivity2 -> Right;
  793.         Connectivity.Top = ChainEdgeLists(Connectivity1 -> Top,
  794.                           Connectivity2 -> Top);
  795.         Connectivity.Bottom = ChainEdgeLists(Connectivity1 -> Bottom,
  796.                          Connectivity2 -> Bottom);
  797.         UpdateBlackHoles(Connectivity1 -> Right, Connectivity2 -> Left,
  798.                  CAGD_CONST_U_DIR);
  799.         break;
  800.     case CAGD_CONST_V_DIR:
  801.         Connectivity.Top = Connectivity1 -> Top;
  802.         Connectivity.Bottom = Connectivity2 -> Bottom;
  803.         Connectivity.Right = ChainEdgeLists(Connectivity1 -> Right,
  804.                         Connectivity2 -> Right);
  805.         Connectivity.Left = ChainEdgeLists(Connectivity1 -> Left,
  806.                            Connectivity2 -> Left);
  807.         UpdateBlackHoles(Connectivity1 -> Bottom, Connectivity2 -> Top,
  808.                  CAGD_CONST_V_DIR);
  809.         break;
  810.     default:
  811.         SYMB_FATAL_ERROR(SYMB_ERR_DIR_NOT_CONST_UV);
  812.     }
  813.     return Connectivity;
  814. }
  815.  
  816. /*****************************************************************************
  817. * DESCRIPTION:                                                               *
  818. * Test for black holes between two adjacent edges and compensate as needed.  *
  819. *                                                                            *
  820. * PARAMETERS:                                                                *
  821. *   L1, L2:  Two edges to search and fill black holes.                       *
  822. *   Dir:     Direction of edge. Either U or V.                               *
  823. *                                                                            *
  824. * RETURN VALUE:                                                              *
  825. *   void                                                                     *
  826. *****************************************************************************/
  827. static void UpdateBlackHoles(EdgeStruct *L1,
  828.                  EdgeStruct *L2,
  829.                  CagdSrfDirType Dir)
  830. {
  831.     while (L1 != NULL && L2 != NULL) {
  832.     int Inc,
  833.         Index1 = L1 -> Index,
  834.         Index2 = L2 -> Index;
  835.     RectangularStruct
  836.         *Rect1 = L1 -> Rect,
  837.         *Rect2 = L2 -> Rect;
  838.  
  839.     Inc = UpdateOneBlackHole(Rect1, Index1, Rect2, Index2, Dir);
  840.  
  841.     if (Inc & 0x01)
  842.         L1 = L1 -> Pnext;
  843.     if (Inc & 0x02)
  844.         L2 = L2 -> Pnext;
  845.     }
  846.  
  847.     if (L1 != NULL || L2 != NULL)
  848.     IritFatalError("Inconsistent edge list in black hole fill");
  849. }
  850.  
  851. /*****************************************************************************
  852. * DESCRIPTION:                                                               *
  853. * Compares the two edges prescribed by the rectangular their are in and the  *
  854. * index within the rectangular for colinearity and close black hole if any.  *
  855. *   Returns two bit flags which edge should be incremented.             *
  856. *                                                                            *
  857. * PARAMETERS:                                                                *
  858. *   Rect1, Index1a:   Index of first edge and rectangular region it is in.   *
  859. *   Rect2, Index2a:   Index of second edge and rectangular region it is in.  *
  860. *   Dir:              Direction of test for black hole. Either U or V.       *
  861. *                                                                            *
  862. * RETURN VALUE:                                                              *
  863. *   int:              Which edge should be incremented. First bit first edge *
  864. *                     second bit second edge.                                *
  865. *****************************************************************************/
  866. static int UpdateOneBlackHole(RectangularStruct *Rect1,
  867.                   int Index1a,
  868.                   RectangularStruct *Rect2,
  869.                   int Index2a,
  870.                   CagdSrfDirType Dir)
  871. {
  872.     int Index1b = (Index1a + 1) % 4,
  873.     Index2b = (Index2a + 1) % 4,
  874.     UVDir = Dir == CAGD_CONST_U_DIR ? 1 : 0;
  875.     CagdRType *P1a, *P1b, *P2a, *P2b, *N1a, *N1b, *N2a, *N2b,
  876.     R1a = Rect1 -> UV[Index1a][UVDir],
  877.     R1b = Rect1 -> UV[Index1b][UVDir],
  878.     R2a = Rect2 -> UV[Index2a][UVDir],
  879.     R2b = Rect2 -> UV[Index2b][UVDir];
  880.  
  881.     if (APX_EQ(R1a, R2a) && APX_EQ(R2a, R2b)) {
  882.     /* Eliminate the trivial case of no black hole. */
  883.     return 0x03;
  884.     }
  885.  
  886.     /* Make sure points are sorted in edge. */
  887.     if (R1a > R1b) {
  888.     SWAP(CagdRType, R1a, R1b);
  889.     SWAP(int, Index1a, Index1b);
  890.     }
  891.     P1a = Rect1 -> Point[Index1a];
  892.     P1b = Rect1 -> Point[Index1b];
  893.     N1a = Rect1 -> Normal[Index1a];
  894.     N1b = Rect1 -> Normal[Index1b];
  895.  
  896.     if (R2a > R2b) {
  897.     SWAP(CagdRType, R2a, R2b);
  898.     SWAP(int, Index2a, Index2b);
  899.     }
  900.     P2a = Rect2 -> Point[Index2a];
  901.     P2b = Rect2 -> Point[Index2b];
  902.     N2a = Rect2 -> Normal[Index2a];
  903.     N2b = Rect2 -> Normal[Index2b];
  904.  
  905.     if (R1a <= R2a && R1b >= R2b) {
  906.     /* Edge 1 contains edge 2. */
  907.     if (!APX_EQ(R2a, R1a)) {
  908.         ProjectPointOnLine(P2a, P1a, P1b, N2a, N1a, N1b);
  909.         MakePointDepend(Rect2, Index2a, Rect1, Index1a, Index1b);
  910.         MakePointDepend(Rect2, Index2a, Rect1, Index1b, Index1a);
  911.         PropagateDependency(Rect2, Index2a);
  912.     }
  913.     if (!APX_EQ(R2b, R1b)) {
  914.         ProjectPointOnLine(P2b, P1a, P1b, N2b, N1a, N1b);
  915.         MakePointDepend(Rect2, Index2b, Rect1, Index1a, Index1b);
  916.         MakePointDepend(Rect2, Index2b, Rect1, Index1b, Index1a);
  917.         PropagateDependency(Rect2, Index2b);
  918.     }
  919.  
  920.     return APX_EQ(R1b, R2b) ? 0x03 : 0x02;
  921.     }
  922.     else if (R2a <= R1a && R2b >= R1b) {
  923.     /* Edge 2 contains edge 1. */
  924.     if (!APX_EQ(R1a, R2a)) {
  925.         ProjectPointOnLine(P1a, P2a, P2b, N1a, N2a, N2b);
  926.         MakePointDepend(Rect1, Index1a, Rect2, Index2b, Index2a);
  927.         MakePointDepend(Rect1, Index1a, Rect2, Index2a, Index2b);
  928.         PropagateDependency(Rect1, Index1a);
  929.     }
  930.     if (!APX_EQ(R1b, R2b)) {
  931.         ProjectPointOnLine(P1b, P2a, P2b, N1b, N2a, N2b);
  932.         MakePointDepend(Rect1, Index1b, Rect2, Index2b, Index2a);
  933.         MakePointDepend(Rect1, Index1b, Rect2, Index2a, Index2b);
  934.         PropagateDependency(Rect1, Index1b);
  935.     }
  936.  
  937.     return APX_EQ(R1b, R2b) ? 0x03 : 0x01;
  938.     }
  939.     else {
  940.     /* Two edges only share part of their domain. Rarely happens. */
  941.     if (GlblFoundPartialShare == FALSE) {
  942.         fprintf(stderr, "Warning: Partial adjacency share detected\n");
  943.         GlblFoundPartialShare = TRUE;
  944.     }
  945.  
  946.     if (R1a < R2a) {
  947.         if (R1b < R2b) {
  948.         ProjectPointOnLine(P2a, P1a, P2b, N2a, N1a, N2b);
  949.         ProjectPointOnLine(P1b, P1a, P2b, N1b, N1a, N2b);
  950.         }
  951.         else {
  952.         ProjectPointOnLine(P2a, P1a, P1b, N2a, N1a, N1b);
  953.         ProjectPointOnLine(P2b, P1a, P1b, N2b, N1a, N1b);
  954.         }
  955.     }
  956.     else {
  957.         if (R1b < R2b) {
  958.         ProjectPointOnLine(P1a, P2a, P2b, N1a, N2a, N2b);
  959.         ProjectPointOnLine(P1b, P2a, P2b, N1b, N2a, N2b);
  960.         }
  961.         else {
  962.         ProjectPointOnLine(P1a, P2a, P1b, N1a, N2a, N1b);
  963.         ProjectPointOnLine(P2b, P2a, P1b, N2b, N2a, N1b);
  964.         }
  965.     }
  966.  
  967.     return R1b > R2b ? 0x02 : 0x01;
  968.     }
  969. }
  970.  
  971. /*****************************************************************************
  972. * DESCRIPTION:                                                               *
  973. * Make Depend point depend on the two Master points.                 *
  974. *                                                                            *
  975. * PARAMETERS:                                                                *
  976. *   DependRect, DependIndex:  Point to become dependent on 1 and 2.         *
  977. *   MasterRect, Master1Index, Master2Index:  The points to depend on.        *
  978. *                                                                            *
  979. * RETURN VALUE:                                                              *
  980. *   void                                                                     *
  981. *****************************************************************************/
  982. static void MakePointDepend(RectangularStruct *DependRect,
  983.                 int DependIndex,
  984.                 RectangularStruct *MasterRect,
  985.                 int Master1Index,
  986.                 int Master2Index)
  987. {
  988.     PointDependStruct
  989.     *PtDepend1 = (PointDependStruct *) IritMalloc(sizeof(PointDependStruct)),
  990.     *PtDepend2 = (PointDependStruct *) IritMalloc(sizeof(PointDependStruct));
  991.  
  992.     PtDepend1 -> Rect = PtDepend2 -> Rect = DependRect;
  993.     PtDepend1 -> Index = PtDepend2 -> Index = DependIndex;
  994.     PtDepend1 -> OtherIndex = Master2Index;
  995.     PtDepend2 -> OtherIndex = Master1Index;
  996.  
  997.     PtDepend1 -> Pnext = MasterRect -> Depend[Master1Index];
  998.     MasterRect -> Depend[Master1Index] = PtDepend1;
  999.     PtDepend2 -> Pnext = MasterRect -> Depend[Master2Index];
  1000.     MasterRect -> Depend[Master2Index] = PtDepend2;
  1001. }
  1002.  
  1003. /*****************************************************************************
  1004. * DESCRIPTION:                                                               *
  1005. * Propagate through dependency information the fact the this point changed.  *
  1006. *                                                                            *
  1007. * PARAMETERS:                                                                *
  1008. *   Rect Index:   "this point".                             *
  1009. *                                                                            *
  1010. * RETURN VALUE:                                                              *
  1011. *   void                                                                     *
  1012. *****************************************************************************/
  1013. static void PropagateDependency(RectangularStruct *Rect, int Index)
  1014. {
  1015.     PointDependStruct *PtDepend;
  1016.  
  1017.     for (PtDepend = Rect -> Depend[Index];
  1018.      PtDepend != NULL;
  1019.      PtDepend = PtDepend -> Pnext) {
  1020.     ProjectPointOnLine(PtDepend -> Rect -> Point[PtDepend -> Index],
  1021.                Rect -> Point[Index],
  1022.                Rect -> Point[PtDepend -> OtherIndex],
  1023.                PtDepend -> Rect -> Normal[PtDepend -> Index],
  1024.                Rect -> Normal[Index],
  1025.                Rect -> Normal[PtDepend -> OtherIndex]);
  1026.     PropagateDependency(PtDepend -> Rect, PtDepend -> Index);
  1027.     }
  1028. }
  1029.  
  1030. /*****************************************************************************
  1031. * DESCRIPTION:                                                               *
  1032. * Update Pt to the closest point on line Pt1Pt2.                 *
  1033. *                                                                            *
  1034. * PARAMETERS:                                                                *
  1035. *   Pt:            Point to be projected on line Pt1Pt2, closing black hole. *
  1036. *   Pt1, Pt2:      The line to project Pt on.                     *
  1037. *   Nrml:          To be computed for Pt from normals at Pt1, Pt2.           *
  1038. *   Nrml1, Nrml2:  Normals at Pt1 and Pt2.                     *
  1039. *                                                                            *
  1040. * RETURN VALUE:                                                              *
  1041. *   void                                                                     *
  1042. *****************************************************************************/
  1043. static void ProjectPointOnLine(CagdPType Pt,
  1044.                    CagdPType Pt1,
  1045.                    CagdPType Pt2,
  1046.                    CagdPType Nrml,
  1047.                    CagdPType Nrml1,
  1048.                    CagdPType Nrml2)
  1049. {
  1050.     CagdRType t;
  1051.     VectorType Vec;
  1052.     PointType NewPt;
  1053.  
  1054.     PT_SUB(Vec, Pt1, Pt2);
  1055.  
  1056.     CGPointFromPointLine(Pt, Pt1, Vec, NewPt);
  1057.     PT_COPY(Pt, NewPt);
  1058.  
  1059.     t = CGDistPointPoint(Pt, Pt2) / CGDistPointPoint(Pt1, Pt2);
  1060.     PT_BLEND(Nrml, Nrml1, Nrml2, t);
  1061.     PT_NORMALIZE(Nrml);
  1062. }
  1063.  
  1064. /*****************************************************************************
  1065. * DESCRIPTION:                                                               *
  1066. * Concatenates two edge lists, destructively.                     *
  1067. *                                                                            *
  1068. * PARAMETERS:                                                                *
  1069. *   L1, L2:   The two edge list to concatenate.                     *
  1070. *                                                                            *
  1071. * RETURN VALUE:                                                              *
  1072. *   EdgeStruct *:   Concatenated edge list.                                  *
  1073. *****************************************************************************/
  1074. static EdgeStruct *ChainEdgeLists(EdgeStruct *L1, EdgeStruct *L2)
  1075. {
  1076.     EdgeStruct
  1077.     *L = L1;
  1078.  
  1079.     if (L1 == NULL)
  1080.     return L2;
  1081.     else {
  1082.     while (L1 -> Pnext != NULL)
  1083.         L1 = L1 -> Pnext;
  1084.     L1 -> Pnext = L2;
  1085.  
  1086.     return L;
  1087.     }
  1088. }
  1089.  
  1090. /*****************************************************************************
  1091. * DESCRIPTION:                                                               *
  1092. * Makes a single triangle from given data.                     *
  1093. * Returns NULL if Triangle is degenerated.                     *
  1094. *                                                                            *
  1095. * PARAMETERS:                                                                *
  1096. *   Pt1, Nrml1, UV1:   Coefficients of first vertex.                 *
  1097. *   Pt2, Nrml2, UV2:   Coefficients of second vertex.                 *
  1098. *   Pt3, Nrml3, UV3:   Coefficients of third vertex.                 *
  1099. *   ComputeUV:         Do we have valid UV information in UV??               *
  1100. *   ComputeNrmls:      Do we have valid normal information in Nrml??         *
  1101. *                                                                            *
  1102. * RETURN VALUE:                                                              *
  1103. *   CagdPolygonStruct *:    Constructed triangle, or NULL if degenerate.     *
  1104. *****************************************************************************/
  1105. static CagdPolygonStruct *MakeTriangle(CagdPType Pt1,
  1106.                        CagdVType Nrml1,
  1107.                        CagdUVType UV1,
  1108.                        CagdPType Pt2,
  1109.                        CagdVType Nrml2,
  1110.                        CagdUVType UV2,
  1111.                        CagdPType Pt3,
  1112.                        CagdVType Nrml3,
  1113.                        CagdUVType UV3,
  1114.                        CagdBType ComputeUV,
  1115.                        CagdBType ComputeNrmls)
  1116. {
  1117.     int i;
  1118.     CagdPolygonStruct *Poly;
  1119.  
  1120.    if (GMColinear3Pts(Pt1, Pt2, Pt3))
  1121.     return NULL;
  1122.  
  1123.     Poly = CagdPolygonNew();
  1124.  
  1125.     PT_COPY(Poly -> Polygon[0].Pt, Pt1);
  1126.     PT_COPY(Poly -> Polygon[1].Pt, Pt2);
  1127.     PT_COPY(Poly -> Polygon[2].Pt, Pt3);
  1128.  
  1129.     if (ComputeUV) {
  1130.     UV_COPY(Poly -> Polygon[0].UV, UV1);
  1131.     UV_COPY(Poly -> Polygon[1].UV, UV2);
  1132.     UV_COPY(Poly -> Polygon[2].UV, UV3);
  1133.     }
  1134.     else {
  1135.     for (i = 0; i < 3; i++)
  1136.         UV_RESET(Poly -> Polygon[i].UV);
  1137.     }
  1138.  
  1139.     if (ComputeNrmls) {
  1140.     PT_COPY(Poly -> Polygon[0].Nrml, Nrml1);
  1141.     PT_COPY(Poly -> Polygon[1].Nrml, Nrml2);
  1142.     PT_COPY(Poly -> Polygon[2].Nrml, Nrml3);
  1143.     }
  1144.     else {
  1145.     for (i = 0; i < 3; i++)
  1146.         PT_RESET(Poly -> Polygon[i].Nrml);
  1147.     }
  1148.  
  1149.     return Poly;
  1150. }
  1151.  
  1152. /*****************************************************************************
  1153. * DESCRIPTION:                                                               *
  1154. * Makes a single edge from given data.                         *
  1155. *                                                                            *
  1156. * PARAMETERS:                                                                *
  1157. *   Rect, Index:    Given the rectangle and the index inside...             *
  1158. *                                                                            *
  1159. * RETURN VALUE:                                                              *
  1160. *   EdgeStruct *:   ...Fetch out the requrested egde.                        *
  1161. *****************************************************************************/
  1162. static EdgeStruct *MakeEdge(RectangularStruct *Rect, int Index)
  1163. {
  1164.     EdgeStruct
  1165.     *Edge = (EdgeStruct *) IritMalloc(sizeof(EdgeStruct));
  1166.  
  1167.     Edge -> Rect = Rect;
  1168.     Edge -> Index = Index;
  1169.     Edge -> Pnext = NULL;
  1170.  
  1171.     return Edge;
  1172. }
  1173.